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Abstract 

We give a cross-disciplinary survey on "population" Monte Carlo al- 
gorithms. In these algorithms, a set of "walkers" or "particles" is used 
as a representation of a high-dimensional vector. The computation is 
carried out by a random walk and split/deletion of these objects. The al- 
gorithms are developed in various fields in physics and statistical sciences 
and called by lots of different terms - "quantum Monte Carlo" , "transfer- 
matrix Monte Carlo", "Monte Carlo filter (particle filter) "/'sequential 
Monte Carlo" and "PERM" etc. Here we discuss them in a coherent 
framework. We also touch on related algorithms - genetic algorithms and 
annealed importance sampling. 



Keywords: quantum Monte Carlo, transfer-matrix Monte Carlo, Monte 
Carlo filter, sequential Monte Carlo, pruned-enriched Rosenbluth method, 
annealed importance sampling, genetic algorithm 
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1 Introduction 



In this paper, we give a cross-disciplinary survey on "population" Monte Carlo 
algorithms. These algorithms, which are developed in various fields, have a com- 
mon structure: A set of "walkers" or "particles" is used for the representation 
of a high-dimensional vector and the computation is carried out by a random 
walk in the state space and split/deletion of these objects. These algorithms 
do not belong to the class of Markov chain Monte Carlo (MCMC, dynamical 
Monte Carlo), although t hey share common features a nd applications. 



Mo nte Ca r lo filter [ Kitagawa |Kitagawa 98| (or sequential Monte 



Carlo jLiu 9S| , |Doucet 01 



Crisan 0C]) algorithm is a topic of recent interest 
in the field of statistical information processing and now being popular in re- 
lated fields, e.g., robot vision [Isard 98]. It is an example of population Monte 
Carlo algorithms defined here. It is, however, not the only example. Population 
Monte Carlo algorithms have proved to be useful tools in a number of fields — 
quantum physics, polymer science, statistical physics, and statistical sciences. 
They are powerful rivals of MCMC in these fields. 

The aim of this paper is neither the design of a novel algorithm nor the pre- 
sentation of new applications. Our goal is to give a minimal set of references and 
explanations for the cross-fertilization among the researchers in different fields. 
It will be useful because few people realize essentially the same algorithms are 



used in both in statistical information processing iKitagawa 96, Kitagawa 98 



Liu 98, Doucet 01 , |Crisan 00 and physics jHctherington 84 , Ccrf 95, Kalos 62 



Ceperley 79, Nightingale 88, Grassberger 0C| |. I also hope that this survey will 



be useful for the development of applications in machine learning and proba- 
bilistic artificial intelligence. 

In Sec. ||, we give a common structure of the population Monte Carlo algo- 
rithms. We also discuss a relation to genetic algorithms and give remarks on 
the origins of the methods. In Sec. [| we give examples of the algorithms both 
in physics and statistics. 



2 An Overview 

2.1 Algorithm — General 

Essentially, the algorithms discussed in this paper are designed for the compu- 
tation of the products of non- negative, sparse, M x M matrices G 1 , G 2 , . . . and 
a non-negative vector X . If we express the vector of nth step as X n then 

X^\i)=Y,G n {i,j)X n (j) (1) 

3 

where 1 < i, j ' < M are indices of components of the vector X n and the matrix 
G™. Here, the dimension M of vector X n is assumed to be very large and we 
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cannot explicitly store the elements of the vector in the memory. Saving the 
storage, we represent the non-negative vector as a weighted superposition of the 
"walker" or "particles" indexed by k (k = 1, . . . , K), each of which is placed at 
j n (k) with a weight w%: 

X n (j)=J2<-S j n {k) , j (2) 

k 

Here, S is the Kronecker delta defined by 

Sm',m =1 (m' = TO) 

= (mV to) (3) 

The algorithm is defined by the iteration of the following steps: 

f . STEP 1: Random walk in the state space. 

Each walker k at the position j n (k) is moved independently to a new 
position i 

j n+1 {k):=i (4) 

according to the probability 

2. STEP 2: Update of the weights. 
For each walker, calculate the factor 

w k =Y,G n a',m) (6) 

i' 

and update the weight of the walker using 

log< +1 :=logW fc +log<. (7) 

3. STEP 3: Reconfiguration. 

• Split walkers with a large weight. Each walker splits into multiple 
walkers, whose total weight is equal to that of the original walker. 

• Remove walkers with a small weight. 

This procedure is called by various different names — reweighting, 
prune/enrichment, reconfiguration, rejuvenation, resampling, 
branching, selection. 

There are several different ways to resample the population of the walkers. 
A simple way is to set 

Q k = =^ (8) 

E/c' w k> 
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and resample the walkers with probability Qk (Note that walkers with 
large Qk can be sampled several times.). Variants are found in the ref- 



erences, for example, Nightingale 88 1. In some variants, the number of 



walkers K is not strictly constant, but fluctuates within a range. 

To complete the algorithm, the initial positions {j°(k)} and weights {u>°} 
of the walkers should be given. In most cases, is a constant independent of 
k and {j°(k)} are independent samples from the initial density defined by X°. 
The choice of X° depends on a specific problem and usually the same as the 
choice of Xq in the corresponding deterministic algorithm based on the direct 
iteration of eq. (fil). For example, we can use an arbitrary density as an initial 



vector X in quantum Monte Carlo (Sec. 3.1), although systematic choice will 
improve the speed of the convergence. In Monte Carlo filtering (Sec. Xq 
represents a prior distribution for initial states of series. 

Let us define [• • •] as an average over random numbers used in the past steps 
of the algorithm. Then, if we assume 

X n (j) = [X n (j)l (9) 

it is not difficult to show that the relation 

X n +\j) = [X n+1 (j)} (10) 

holds with STEP 1 ~ STEP3. This property is essential for the adequacy of 
population Monte Carlo algorithm. When we modify the STEP 3, we should be 
careful to conserve this property. With the condition eq. (|l(]), we expect the 
convergence of the algorithm in the limit of the number of the walkers K — > oo 
with a fixed number n of the iteration []. We refer to [ Crisan OOf for a recent 



result on the convergence of population Monte Carlo algorithms n. Practically, 
the rate of the convergence severely depends on the problem treated by the 
method. 

Finally, we note that it is easy to generalize the algorithm to the cases with 
continuous state space, where the matrix G is replaced by an operator and the 
summations become integrals. Such continuous versions of the algorithm are 
successfully used in some of the examples discussed below, e.g., applications to 
quantum many-body problems. 

2.2 The Role of STEP 3 

The STEP 3 of the algorithm is not necessary for the validity of the algorithm, 
i.e., the condition that we discussed in the previous section is satisfied only with 

1 Note that the convergence is not assured in the limit n — > oo with a fixed number K of 
walkers. It is an imp ortant remark for the problems where the limit n — -> oo is required (for 



example, see Sec. 3.1). 
2 Some rem 



in [ Hctherington 84 , Cerf 95 1 



bv ph ysicists on the property of the algorithm are found 
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STEP 1 and STEP 2. The role of STEP 3 is to suppress the variation of weights 
{wk} and improve the efficiency of the algorithm. 

For some problems, the algorithm without STEP 3 is sufficient. Hereafter, 
we call the algorithm without STEP 3 as "simple" population Monte Carlo. 
Note that without STEP 3, there is no interaction between walkers. Then, 
parallel evolution of the walkers can be replaced by sequential runs each of 
which corresponds to the simulation of a walker. 

In some references, a p art of STEP 3 is separately discussed as a "population 
control" procedure (e.g., |Kalos 62 , Grassbergcr 98~| ) and/or a part of STEP 3 
is merged into STEP 1 (e.g., [ Kalos 62 1). In such cases, we should be careful to 
judge whether the algorithm effectively contains STEP 3 or not. 



2.3 Relation to Genetic Algorithm 

There is an obvious analogy between population Monte Carlo and genetic al- 
gorithms (GA) [ Higuchi "97| . However, there is an essential difference in the 
goal of the algorithm. Population Monte Carlo is a tool for the computation 
of the product of matrices, multiple summations (integrals) and calculation of 
marginals. On the other hand, GA is a tool specialized in optimization. Usually, 
we can modify population Monte Carlo for the optimization. However, the con- 
verse is often not true. Specifically, we should be careful to the introduction of 
crossover operators and other tricks popular in GA to population Monte Carlo, 
because they can easily spoil the convergence to the exact result in the limit of 
infinite number of walkers, K —> oo. 

There are some theoretical studies on GA that assume the absence of the 
crossover operation. It might be interesting to apply them to the study of 
population Monte Carlo. 



2.4 Origins 

It is difficult to fix the origin of the population Monte Carlo. An origin seems the 
algorithms for the calculat ion of the elements of an inver se Matrix developed by 
von Neumann and others [ Hammersley 64 , Drcitlcin 85||, wh ich lacks the STEP 
3. Another reference is Metropolis and Ulam [ Metropolis 49 fj], which discuss a 
solver of the Schrodinger equation by simulation. The algorithm is attributed to 
Fermi (It is an origin of quantum Monte Carlo algorithms with random walkers 
discussed in Sec. 3.1). The algorithm for eigenvalue problems are also discussed 
in the classic book Hammersley 64 on Monte Carlo. Simple algorithm for self- 
avoiding walks is introduced by Rosenbluth and Rosenbluth [Rosenbluth 55 1. 
According to [Crisan 00 1, one of the earliest studies of simple algorithm in sta- 
tistical science is found in [Handschin 69 1. 



3 It is interesting to point out that this paper by famous physicists appeared in the Journal 
of American Statistical Association, which is one of the well-known journals of statistics. 
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It is also not easy to determine who introduced STEP 3. For example, 
STEP 3 is included in the algorithm in Metropolis and Ulam [Metropolis 49 1 
and Kalos [Kalos 62 in a somewhat implicit manner. The "enrichment" al- 
gorithm [ Wall 59{ for self-avoiding walks contains a very special case of STEP 
3. Anyway, many authors have used the algorithms with STEP 3 by the end of 
1980s as a version of "quantum Monte Carlo" or "transfer-matrix Monte Carlo" , 
while the algorithms without STEP 3 have also been used up to now. 

We left a comprehensive treatment on the history of the population Monte 
Carlo methods to future surveys. 



3 Examples 



3.1 Quantum Many-Body Problems 



First, we discuss "qu a ntum Monte Carlo" algorithms [ ! | Ceperley 79 



Schmidt 84| jSchmidt 92| |De R.aedt 92j, JHctherington 84j |Cerf 95| , |Negele 86 
Anderson 7q , Blankcnbeclcr 83 , John 87 1. In spite of their name, they are not 
algorithms specialized to quantum mechanics. They are most naturally under- 
stood as methods for the approximation of the smallest (or largest) eigenvalue 
and the corresponding eigenvector of a large sparse matrix. Essentially, they 
are stochastic versions of the power method for eigenvalue problems. 

Consider a symmetric, non-negative matrix A. To formulate the eigenvalue 
problem into population Monte Carlo, we introduce the exponential exp(— f3A) 
of A. For sufficiently large /?, the expression 



exp(-/L4) X° 



(11) 



approximately gives an eigenvector with the smallest eigenvalue for an arbitrary 
vector X° that has non-zero projection to the eigenvector. In general, the 
calculation of exp(— f3A) is difficult. In most of interesting problems, however, 
A is decomposed into several simple components, e.g., A = B + C where B and 
C is symmetric, non-negative matrices whose exponential is easily calculated. 
Thus, the exponential of A is expressed as 



exp(— (3 A) = lim 



N- 



ex P(^^-- B )exp(-^:C) 



A' 



(12) 



It is straightforward to fit the problem into the form of the equation eq. (Ill) with 



G n = exp(-|s)exp(-Ac) 



(13) 



4 There is a different type of "quantum Monte Carlo" algorithms which are based on MCMC 
(dynamical Monte Carlo). They simulate a "world-line" instead of a collection of random 
walkers. 
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where N is a sufficiently large number and n = 0, . . . , N — 1. 

Another approach to the eigenvalue problem is a method based on an itera- 
tive solver of the eigenvalue equation. Consider the product of matrices 



G n =B- 1 (XI -D) 



(14) 



where A = B + D and D are a diagonal matrix (7 denotes an identity matrix.). 
Under the suitable condition on the "Green's function" B^ 1 and a properly 
chosen value of A 0, the vector X n converges to an eigenvector of A. 

In the treatment on quantum mechanical problems, the matrix A is usually 
Hamiltonian of the given system. It is easy to construct continuous version of 
the algorithm, where the matrices and summation is replaced by operators and 
integrals, respectively. In such cases, the operator B is often represented by the 
summation of Laplace operators and a constant. 

These algorithms are successfully applied to various problems in quantum 



physics, e.g., the computation of the properties of nuclei |Kalos 62, Negele 



supcrfluids [Ceperley 79, Schmidt 92 1, quantum spin systems [ Blankenbecler 83 



De Raedt 92 , and quantum dynamics [John 87 1. They are called by the terms 
"projector Monte Carlo" , "Green's function Monte Carlo" , and "diffusion Monte 
Carlo" |. 

In some cases, we should deal with matrices or eigenvectors that are not non- 
negative. Formally, these cases can be dealt with the introduction of negative 
weight of walkers. It, however, seriously spoils the efficiency and convergence 
of the algorithm. It is an example of "negative sign crisis" in stochastic compu- 
tation of quantum mechanical problems. Although there are some interesting 



ideas to improve the algorithm [Kalos 84, Carlson 85, Negele 89, Anderson 91 
Liu 94 1, it seems difficult to remove the difficulty entirely with a smart trick, 



because it originates from radical difference between classical stochastic systems 
and quantum systems. 



3.2 Lattice Spin Systems 

Using the transfer matrix (transfer integral) formalism, we can translate a calcu- 
lation for classical statistical mechanics (and combinatorics) to the computation 
of a product of matrices. Then, it is a natural idea to apply population Monte 
Carlo algorithms to these problems. For example, consider a classical spin model 
(Markov field model) defined on a L\ x L2 (L2 >> L\) strip We define a set 
of subsystems each of which consists of the 1 , . . . , n rows of the original strip 

5 When we implement it in a form of population Monte Carlo, the choice of the constant A is 
absorbed in STEP 2 and STEP 3. Then we can obtain the eigenvalue A and the corresponding 
eigenvector simultaneously. 

6 Here we avoid the details of the terminology, because the classification and naming of 
these algorithms may depend on the authors. 

7 For a three-dimensional problem, a L\ X L\ X L2 column substitutes for the strip. 
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(n = 1,2, . . . , L2). Hereafter we assume nearest neighbor interaction on the lat- 
tice, although extensions to cases with next-nearest interaction etc. are formally 
easy. The partition function Z n (S n ) of the subsystem n conditioned with the 
values of the variables in the nth 

EE"" E exp(-f3E n (S 1 ,S 2 ,...,S n )) (15) 

Si S2 S n -i 

where f3 is inverse temperature and E n (Si, S2, • ■ ■ , S n ) is the energy of the nth 
subsystem. When each variable Si n is a binary variable, the number of the 
possible values of S n is 2 Ll . Then, a recursion relation 

Z n +i{S n +i) = 

E s „ Z n (S n ) exp(-(3AE n (S n , S n+1 )) (16) 

where 

AE n (S n , S n +i) = 
E n +i(Si, . . . , S n , S n +i) — E„(Si, . . . , S n ) (17) 

holds. Obviously, it is essentially the same as the filtering formula for general 
state space models and the recursive formula of the Baum- Welch algorithm for 
hidden Markov models. We define Z a = 1 and AE (Si) = Ei(S\) for the free 
boundary condition at the top of the strip (Here and hereafter, the symbol So 
in the formulae should be neglected.). 

A population Monte Carlo algorithm is constructed by setting G and X as 

G n (S n+1 , S n ) = exp(~(3AE n (S n , S n+1 )) (18) 

X n (S n ) = Z n (S n ). (19) 

The algorithm is called "tran sfer-matrix Mo n te Carlo" an d applied to classical 
Heisenberg and XY models 
the number of Penrose tiling 



Nightingale 88j , Thijsscn 9C ] and enumeration of 



Shaw 91 on a lattice. 



3.3 Bayesian Computation: MCF and Sequential Monte 
Carlo 

Algorithms that are similar to the "transfer-matrix Monte Carlo" are developed 
by statisticians for the models of time-series analysis and models with sequential 
structures (e.g., models for gene-propagation analysis). In these algorithms, we 
consider marginal likelihood conditioned with the state x n of the latest step 
(time-step) n 

z„( Xn ) - EE-- - E H*o)p(xo\xi) 

Xq X\ X n — \ 

P{yi\x\) ■ ■ ■ P(x n - 1 \x n )P(y n \x n ) 

(20) 
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instead of partition function []. Here P{y\x) is likelihood associated with data 
y (representation of observational noise). P(x\x') and P(x) are prior probabil- 
ity for the state transition x — > x' (representation of system noise) and prior 
probability of initial state x, respectively. The application of population Monte 

Carlo algorith ms is straightf o rward. 

Kitagawa [ Kitagawa 96| , Kitagawa 98 1 introduced "Monte Carlo filter" 
(MCF, "particle filter") for non-Gaussian/non-linear state space models. Gor- 
don and coworkers independently developed an idea similar to the one by Kita- 
gawa, which they call "bootstrap filter" | Gordon 95 , Gordon 95| . Another con- 
tribution is "sequential imputation" that first used in gene-propagation analy- 
sis [ Irwin 94 1 and later used for blind deconvolution [ Liu 95f and other problems. 
The original form of sequential imputation is a simple algorithm that lacks STEP 
3, but STEP 3 is introduced later ("rejuvenated sequential imputation"). The 
term "sequential Monte Carlo" is also used to indicate a f amily o f algorithm s 
that include the above-mentioned methods as special cases [Liu 98, Crisan 00 1. 
In [Liu 98, Crisan 0C], more references in statistical sciences, some of which 
are earlier than above-mentioned studies, are discussed. Applications in robot 
vision are described in [[sard 98]. A forthcoming book [Doucet 01 on sequen- 



tial Monte Carlo by statisticians will be a good guide to the applications in 
statistical sciences. 



3.4 Polymer Science 

"Polymer" means a flexible chain consists of small molecules ("monomers") 
connected in a fixed order. It is called a heteropolymer (<-> homopolymer) 
when it is made by monomers of more than two species. Well-known examples 
of heteropolymer is protein, RNA and DNA. A good model of polymer is a 
flexible chain where each unit is mutually repulsed by short range force (van 
der Waals core) and repulsed/attracted by longer range force. In lattice models, 
the former interaction is modeled by "self-avoiding condition" , i.e., the condition 
that a lattice point should not be occupied by more than two monomers. 

Population Monte Carlo algorithms for the study of polymer models are 
designed in a manner similar to those for lattice spin models. Subchains of length 
n = 1, 2, . . . , N that consists of the first n monomers of the original polymer of 
length N is used as subsystems for the definition of the algorithm. An important 
difference is that the index of a vector X n is not the coordinates r„ of the 
monomer located at the end of a subchain, but coordinates {r' n }, [n! = 1, . . . , n) 
of all monomers in the subchains. This is because a pair of monomers with a 
large distance along the chain can interact each other, even with short range 
interaction. Thus, the recursion is better described by the Gibbs distribution 
P n of subchains, rather than partition functions Z n . 

8 We can also use formulati on based on posterior distributions of subsystems, which is 
similar to the one used in Sec. 3.4 It is useful in some problems, e.g., smoothing of time- 
series. 



If we express the Gibbs distribution of a subchain n as P n (ri, r%, . . . , r n ), 
the recursion relation 



P n +i(n,r 2 , . . . ,r n+ i) oc 
P n {ri,r 2 , . ■ . ,r n ) ■ exp(-f3AE({rx, 



, r n}, r n+l)) 



(21) 



holds, where Ai?({ri, . . . , r n }, r n +i) is the interaction energy between a new 
monomer and the monomers in the subchain n. For lattice models, we set 
AE = +00 if the conformation of the subchain specified by {ri, . . . , r„+i} is 
non-physical, i.e., self-avoiding condition or connectivity condition of the chain 
is violated. 

To define a population Monte Carlo algorithm, we must specify the index i 
of X n (i) in the equation eq. (|l|). The index i should determine the conformation 
{r\, . . . , r l n } of the nth subchain and should have the dimension N independent 
of n. Here we define i as a vector with components i a = r*(s < n) and i a = 
l(s > n), i.e., 



15 ' 25 



1,1, 



,1) 



(22) 



where 1 is the coordinates of an arbitrary point, for example, the origin. Next 
we define X n (i) as 

X n {i)<xP n (r{,ri,...,ri) (23) 

for i of the form eq. ( |22| ) and set X n (i) = for other indices i. Finally, we 
define the matrix G n as 



G" =exp(-/3A£({rt,. 



3s) III 



n+ 2^s - 1). 



(24) 



With these expressions, we can write the recursion eq. ( pi| ) in the form of the 
equation eq. ([j]), which gives population Monte Carlo algorithms based on a 
simulation of growth/selection of polymer subchains (here we start from X 1 
instead of X°, with which the index n equals to the length of the subchain). 
Population Mont e Carlo algorithms f or po lymer models are developed by 



two or more groups ([ Garel 90 , |Orland 98 ] and [ Grassberger 97 , Grassbcrger 98 



Bastolla 9g| , Grassbcrger 00 1 ) [J Grassberger pointed out that his algorithm, 



PERM (pruned-enriched Rosenbluth method) is understood as an unification 
of "enri chment algorithm" |Wall 5E ] (STEP 1 and 3) and "Rosenbluth algo- 
rithm" [|Rosenblutn~55| , |Seno 96|| (STEP 1 and 2). PERM has proved to be an 



useful tool to study statistical mecha nics of long homopolymers jGrassbergcr 97 
and other lattice statistical models [Grassberger 98, Grassberger 0C|. It is also 
the most efficient algorithm to compute finite temperature properties of lattice 
protein models Bastolla 98, Grassberger 00 1, although it is recently challenged 
by a MCMC algorithm developed by the author and coworkers | Chikcnji 99[ . 



9 A related problem (counting 
Monte Carlo. 



of "meanders") is treated in [Golinelli 00] by population 



10 



3.5 Bayesian Computation: Annealed Importance Sam- 
pling 



In the example of Sec. 3.1, an intermediate state indexed by n is a fictitious 



one and we are only interested in the limit of large n. In Sec. 3.2-3.4, we dealt 
with examples where the index n corresponds to an index in the real world, 
e.g., the index of row of a lattice, time-step of series, and the length of the 
subchain. Here we consider cases where the index n indicates a system with a 
parameter j n , e.g., temperature in statistical physics and hyperparameters in 
Bayesian models. Such an example is discussed by Neal and called by "annealed 



importance sampling" [Neal 98 1 



Consider a parametric family of distributions {P 7n (x), n — 1, .. . , N}. An 
example is 

P yn (x)={P(y\x)}^P(x) (25) 

where P(y\x) is likelihood with data y and P(x) is a prior. We can formally 
define a "recursion relation" by 

p ^ x) = ^M- p -'^ x) - (26) 

Unfortunately, population Monte Carlo algorithm derived from eq. ( p6| ) is trivial 
and inefficient, because the matrix G n (x, x') defined by the recursion is diagonal. 

An idea by Neal is to introduce MCMC steps into the algorithm. The inser- 
tion of probabilistic state-changes that makes P 7ti (x) invariant corresponds to 
the insertion of a non-diagonal matrix G n to the right-hand-side of the equa- 
tion cq. (|l|), which satisfies 

X n (i) = J2G n (i,j)X n (j). (27) 



Obviously, the introduction of such G n does not have an effect on the validity 
of the results. The adequacy of the procedure is also ensured by the relation 

X n (j) = [X n (j)) (28) 

is not changed by the addition of MCMC steps. The resultant algorithm consists 
of finite MCMC sweeps at each n and STEP 2 with 

where x = j(k) is the position of the fcth walker. 

As far as I know, Neal does not discuss algorithms with STEP 3, but there 
seems no reason against the use of STEP 3 in the algorithm. On the other 
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hand, as Neal himself pointed out, the algorithm without STEP 3 can be re- 
garded as a version of non-adiabatic thermodynamic integration discussed by 
Jarzynski [Jarzynski 97, Ecndrix 01 1. 



It is interesting to remark that the introduction of MCMC steps [Crisan 00 1 
is possible for some of the other population Monte Carlo algorithms. 
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